package scu.maqiang.elastodynamics;

import java.io.FileWriter;
import java.io.IOException;
import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES1L21;
import scu.maqiang.fes.FES1L2PE;
import scu.maqiang.homogenization.SOTS1DL2PE;
import scu.maqiang.mesh.MatrixFunc;
import scu.maqiang.mesh.Mesh1L2;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.DirectSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

public class AxisPEDynamicsSOTS1DL2 {

	public static void main(String[] args) {
		long begin = System.nanoTime();
		int N = 32;
		double ep = 1.0 / N;
		
		ScalarFunc coefG1 = (xyz, label, param) -> label == 10 ? 20 : 4;
		ScalarFunc coefG2 = (xyz, label, param) -> label == 10 ? 30 : 6;
		ScalarFunc coefLambda = (xyz, label, param) -> label == 10 ? 15 : 3;
		ScalarFunc coefEp1 = (xyz, label, param) -> label == 10 ? 1.0e-5 : 0.02e-5;
		ScalarFunc coefE1 = (xyz, label, param) -> label == 10 ? 5.0 : 0.3;
		ScalarFunc coefE2 = (xyz, label, param) -> label == 10 ? 6.0 : 0.4;
		ScalarFunc coefrho = (xyz, label, param) -> label == 10 ? 1 : 0.5;
		ScalarFunc coeff1 = (xyz, label, param) -> label == 10 ? 1 : 1;
		ScalarFunc coefrf = (xyz, label, param) -> label == 10 ? 1 : 1;

		double[][] D1 = {{20, 15, 5.0}, {15, 30, 6.0}, {-5.0, -6.0, 1.0e-5}};
		double[][] D2 = {{4, 3, 0.3}, {3.0, 6.0, 0.4}, {-0.3, -0.4, 0.02e-5}};
		MatrixFunc matFunc = (xyz, label, param) -> label == 10? D1: D2;
//		ScalarCoefFunc coefG1 = (xyz, label, param) -> label == 10 ? 20 : 4;
//		ScalarCoefFunc coefG2 = (xyz, label, param) -> label == 10 ? 30 : 6;
//		ScalarCoefFunc coefLambda = (xyz, label, param) -> label == 10 ? 15 : 3;
//		ScalarCoefFunc coefEp1 = (xyz, label, param) -> label == 10 ? 1.0e-5 : 0.02e-5;
//		ScalarCoefFunc coefE1 = (xyz, label, param) -> label == 10 ? 5.0e-2 : 0.3e-2;
//		ScalarCoefFunc coefE2 = (xyz, label, param) -> label == 10 ? 6.0e-2 : 0.4e-2;
//		ScalarCoefFunc coefrho = (xyz, label, param) -> label == 10 ? 1 : 0.5;
//		ScalarCoefFunc coeff1 = (xyz, label, param) -> label == 10 ? 1 : 1;
//		ScalarCoefFunc coefrf = (xyz, label, param) -> label == 10 ? 1 : 1;
		//���㵥������
		Mesh1L2 cellMesh = new Mesh1L2();
		cellMesh.line(0, 1.0, 256);
		ScalarFunc cellRegionFunc1D = (x, llable, param) -> {
			return Math.abs(x[0] - 0.5) < 0.25?1.0:-1.0;
		};
		cellMesh.setDomainLabel(cellRegionFunc1D, null, 10);
		
		FES1L2PE fsCell = new FES1L2PE(cellMesh);
		double[][] FOCS = new double[3][fsCell.GetNdof()];
		double[][] SOCS = new double[6][fsCell.GetNdof()];
		double[][] D = new double[3][3];
		SOTS1DL2PE ss = new SOTS1DL2PE(fsCell);
		ss.AxisPECellSolution(new ScalarFunc[] { coefG1, coefG2, coefLambda, 
				                                     coefEp1, coefE1, coefE2, 
				                                     coefrho, coeff1, coefrf }, null, FOCS, D, SOCS);
		System.out.println("OK");
		
		System.out.println(MVO.toString(D));
		
		ss.computeHomogenizedAxisPE(new double[] {20, 30.0, 15.0, 1.0e-5, 5.0, 6.0}, 
				                    new double[] {4, 6, 3, 0.02e-5, 0.3, 0.4}, 0.5);
		double time = 0.0, dt = 0.005, endTime = 1.0;
		double alpha = 0.25, delta = 0.5;
		double c0 = 1.0 / (alpha * dt * dt);
		double c1 = delta / (alpha * dt);
		double c2 = 1.00 / (alpha * dt);
		double c3 = 0.50 / alpha - 1.0;
		double c4 = delta / alpha - 1.0;
		double c5 = 0.5 * dt * (delta / alpha - 2.0);
		double c6 = dt * (1 - delta);
		double c7 = delta * dt;

		Mesh1L2 homoMesh = new Mesh1L2();
		homoMesh.line(0, 1, 4096);
		homoMesh.translate(1.0);
		FES1L2PE fsHomo = new FES1L2PE(homoMesh);
		SRMatrix KHomo = new SRMatrix(fsHomo.GetNdof());
		//SRMatrix M0Homo = new SRMatrix(fsHomo.GetNdof());
		SRMatrix MHomo = new SRMatrix(fsHomo.GetNdof());
		double[] RHSHomo = new double[fsHomo.GetNdof()];
		double[] uHomo = new double[fsHomo.GetNdof()];
		double[] vHomo = uHomo.clone();
		double[] aHomo = uHomo.clone();
		double[] u0Homo = uHomo.clone();
		double[] v0Homo = uHomo.clone();
		double[] a0Homo = uHomo.clone();
		
		FES1L21 fsHomo1 = new FES1L21(homoMesh);
		SRMatrix M0Homo = new SRMatrix(fsHomo1.GetNdof());
		double[] RHSHomo1 = new double[fsHomo1.GetNdof()];
		
		double[][] DHomo = {{D[0][0], D[0][2], D[1][1]}, 
				            {D[0][2], D[0][1], D[1][2]}, 
				            {-D[1][1], -D[1][2], D[1][0]}};
		
		fsHomo.assemblePEStiff(new double[] { D[0][0], D[0][1], D[0][2], D[1][0], D[1][1], D[1][2]}, BVPType.AXISYMMETRIC, KHomo);
		//fsHomo.assemblePressure(new double[] {-1.0}, BVPType.AXISYMMETRIC, RHSHomo, 2);
		//KHomo.sort();
		//System.out.println(KHomo);
//		fsHomo.assemblePEStiff(new double[] { D[0][0], D[0][1], D[0][2] }, BVPType.COMMON, KHomo);
		//fsHomo.assemblePEMass(new double[] {1.0}, BVPType.AXIS_CONSIST_MATRIX, M0Homo);
		fsHomo.assemblePEMass(new double[] {D[2][0]}, BVPType.AXIS_CONSIST_MASS, MHomo);
		
		SRMatrix LHomo = KHomo.add(c0, MHomo);
		double[] uvaHomo = uHomo.clone();
		
		//ScalarCoefFunc initFunc =  (x, label, param) -> 0.01 * Math.sin(0.5 * Math.PI * (x[0] - 1.0)); 
		ScalarFunc initFunc =  (x, label, param) -> 0.001 * (x[0] - 1.0); 
		//ScalarCoefFunc initFunc =  (x, label, param) -> 0.0; 
		double[] uu0Homo = fsHomo.valueFromFunc(initFunc, null);

		//System.out.println(Arrays.toString(uu0Homo));
		for(int i = 0; i < uu0Homo.length; i++) {
			u0Homo[2 * i] = uu0Homo[i];
		}
		double[] Ku0Homo = KHomo.mul(u0Homo);
		
		fsHomo1.assembleMass(new double[]{D[2][0]}, BVPType.AXIS_CONSIST_MASS, M0Homo);
		//fsHomo1.assembleHeatFlux(new double[] {-1.0e2}, BVPType.AXISYMMETRIC, RHSHomo1, 2);
		for(int i = 0; i < fsHomo1.GetNdof(); i++) {
			RHSHomo1[i] = RHSHomo1[i] - Ku0Homo[2 * i];
		}
		
		DirectSSolver mHomoSolver = new DirectSSolver(M0Homo);
		double[] a0Homo1 = new double[fsHomo1.GetNdof()];
		a0Homo1 = mHomoSolver.CholeskySolve(RHSHomo1, 1);
		for(int i = 0; i < fsHomo1.GetNdof(); i++) {
			a0Homo[2 * i] = a0Homo1[i];
		}
		Arrays.setAll(aHomo, i -> a0Homo[i]);
		
		
		for(int i = 0; i < fsHomo.GetNdof(); i++) {
			uvaHomo[i] = c0 * u0Homo[i] + c2 * v0Homo[i] + c3 * a0Homo[i];
		}
		MHomo.mul(uvaHomo, RHSHomo);
//		String outputFile1 = "PEFreeVibrationSOTS1Du.dat";
//		String outputFile2 = "PEFreeVibrationSOTS1Dv.dat";
//		time = 0.0;
//		
//		double[][] uvtimeHomo = new double[2][homoMesh.getNv()];
//		fsHomo.extractComponent(u0Homo, uvtimeHomo);
//		homoMesh.toTecplot(outputFile1, time, uvtimeHomo[0]);
//		homoMesh.toTecplot(outputFile2, time, uvtimeHomo[1]);
		
		//fsHomo.applyBC_MBN(LHomo, Direct.X, 1);
		fsHomo.applyBC_MBN(LHomo, Direct.Y, 1, 2);
		
		DirectSSolver homoSolver = new DirectSSolver(LHomo);


		Mesh1L2 fineMesh = new Mesh1L2();
		fineMesh.line(0, 1, 4096);
		fineMesh.translate(1.0);
		ScalarFunc regionFineFunc1D = (x, llable, param) -> {
			double xCell = x[0] * N - (int) (x[0] * N);
			return Math.abs(xCell - 0.5) < 0.25 ? 1.0 : -1.0;
		};
		fineMesh.setDomainLabel(regionFineFunc1D, null, 10);
		
		FES1L2PE fsFine = new FES1L2PE(fineMesh);
		SRMatrix KFine = new SRMatrix(fsFine.GetNdof());
		SRMatrix MFine = new SRMatrix(fsFine.GetNdof());
		double[] RHSFine = new double[fsFine.GetNdof()];
		double[] uFine = new double[fsFine.GetNdof()];
		double[] vFine = uFine.clone();
		double[] aFine = uFine.clone();
		double[] u0Fine = uFine.clone();
		double[] v0Fine = uFine.clone();
		double[] a0Fine = uFine.clone();
		
		FES1L21 fsFine1 = new FES1L21(fineMesh);
		SRMatrix M0Fine = new SRMatrix(fsFine1.GetNdof());
		double[] RHSFine1 = new double[fsFine1.GetNdof()];
		
		fsFine.assemblePEStiff(new ScalarFunc[] { coefG1, coefG2, coefLambda, coefEp1, coefE1, coefE2 }, null,
				BVPType.AXISYMMETRIC, KFine);
		
		//KFine.sort();
		//System.out.println(KFine);
		fsFine.assemblePEMass(new ScalarFunc[] { coefrho }, null, BVPType.AXIS_CONSIST_MASS, MFine);
		SRMatrix LFine = KFine.add(c0, MFine);
		double[] uvaFine = uFine.clone();
		
		double[] uu0Fine = fsFine.valueFromFunc(initFunc, null);
		
		for(int i = 0; i < uu0Fine.length; i++) {
			u0Fine[2 * i] = uu0Fine[i];
		}
		double[] Ku0Fine = KFine.mul(u0Fine);
		
		fsFine1.assembleMass(new ScalarFunc[]{coefrho}, null, BVPType.AXIS_CONSIST_MASS, M0Fine);
		//fsFine1.assembleHeatFlux(new double[] {-1.0e2}, BVPType.AXISYMMETRIC, RHSFine1, 2);
		for(int i = 0; i < fsFine1.GetNdof(); i++) {
			RHSFine1[i] = RHSFine1[i] - Ku0Fine[2 * i];
		}
		
		DirectSSolver mFineSolver = new DirectSSolver(M0Fine);
		double[] a0Fine1 = new double[fsFine1.GetNdof()];
		mFineSolver.LUSolve(RHSFine1, 1, 0.04);
		for(int i = 0; i < fsFine1.GetNdof(); i++) {
			a0Fine[2 * i] = a0Fine1[i];
		}
		Arrays.setAll(aFine, i -> a0Fine[i]);
		
		//for(int i = 0; i < fsFine.GetNdof(); i++) {
		//	System.out.println(RHSFine[i]);
		//}
		//System.out.println(Arrays.toString(RHSHomo));
		//System.out.println(Arrays.toString(RHSFine));
		
		//System.out.println(Arrays.toString(a0Fine));
		//System.exit(0);
		
		for(int i = 0; i < fsFine.GetNdof(); i++) {
			uvaFine[i] = c0 * u0Fine[i] + c2 * v0Fine[i] + c3 * a0Fine[i];
		}
		MFine.mul(uvaFine, RHSFine);
		time = 0.0;
		//fsFine.applyBC_MBN(LFine, Direct.X, 1);
		fsFine.applyBC_MBN(LFine, Direct.Y, 1, 2);
		
		DirectSSolver fineSolver = new DirectSSolver(LFine);

		double[] fineU0 = new double[fsFine.GetNdof()];
		double[] fineU1 = new double[fsFine.GetNdof()];
		double[] fineU2 = new double[fsFine.GetNdof()];
		
		double[][] fineU_uphi = new double[2][fineMesh.getNv()];
		double[][] fineU0_uphi = new double[2][fineMesh.getNv()];
		double[][] fineU1_uphi = new double[2][fineMesh.getNv()];
		double[][] fineU2_uphi = new double[2][fineMesh.getNv()];
		double[][] u012 = new double[4][];
		double[][] phi012 = new double[4][];
		
		String outputFile1f = "AxisPED\\SOTS1Duf_"+N+".dat";
		String outputFile10 = "AxisPED\\SOTS1Du0_"+N+".dat";
		String outputFile11 = "AxisPED\\SOTS1Du1_"+N+".dat";
		String outputFile12 = "AxisPED\\SOTS1Du2_"+N+".dat";
		String outputFile2f = "AxisPED\\SOTS1Dphif_"+N+".dat";
		String outputFile20 = "AxisPED\\SOTS1Dphi0_"+N+".dat";
		String outputFile21 = "AxisPED\\SOTS1Dphi1_"+N+".dat";
		String outputFile22 = "AxisPED\\SOTS1Dphi2_"+N+".dat";
		
		String outputFileD1f = "AxisPED\\SOTS1DuDf_"+N+".dat";
		String outputFileD10 = "AxisPED\\SOTS1DuD0_"+N+".dat";
		String outputFileD11 = "AxisPED\\SOTS1DuD1_"+N+".dat";
		String outputFileD12 = "AxisPED\\SOTS1DuD2_"+N+".dat";

		String outputFileD2f = "AxisPED\\SOTS1DuDphif_"+N+".dat";
		String outputFileD20 = "AxisPED\\SOTS1DuDphi0_"+N+".dat";
		String outputFileD21 = "AxisPED\\SOTS1DuDphi1_"+N+".dat";
		String outputFileD22 = "AxisPED\\SOTS1DuDphi2_"+N+".dat";

		String outputFileDDf = "AxisPED\\SOTSDDf_"+N+".dat";
		String outputFileDD0 = "AxisPED\\SOTSDD0_"+N+".dat";
		String outputFileDD1 = "AxisPED\\SOTSDD1_"+N+".dat";
		String outputFileDD2 = "AxisPED\\SOTSDD2_"+N+".dat";
		
		
		fsFine.extractComponent(uFine, fineU_uphi);
		u012[0] = fineU_uphi[0];
		u012[1] = fineU_uphi[0];
		u012[2] = fineU_uphi[0];
		u012[3] = fineU_uphi[0];
		
		phi012[0] = fineU_uphi[1];
		phi012[1] = fineU_uphi[1];
		phi012[2] = fineU_uphi[1];
		phi012[3] = fineU_uphi[1];
		fineMesh.toTecplot(outputFile1f, time, fineU_uphi[0]);
		fineMesh.toTecplot(outputFile10, time, fineU_uphi[0]);
		fineMesh.toTecplot(outputFile11, time, fineU_uphi[0]);
		fineMesh.toTecplot(outputFile12, time, fineU_uphi[0]);
		fineMesh.toTecplot(outputFile2f, time, fineU_uphi[0]);
		fineMesh.toTecplot(outputFile20, time, fineU_uphi[0]);
		fineMesh.toTecplot(outputFile21, time, fineU_uphi[0]);
		fineMesh.toTecplot(outputFile22, time, fineU_uphi[0]);
		fineMesh.toTecplot(outputFileD1f, time, fineU_uphi[0]);
		fineMesh.toTecplot(outputFileD10, time, fineU_uphi[0]);
		fineMesh.toTecplot(outputFileD11, time, fineU_uphi[0]);
		fineMesh.toTecplot(outputFileD12, time, fineU_uphi[0]);
		fineMesh.toTecplot(outputFileD2f, time, fineU_uphi[0]);
		fineMesh.toTecplot(outputFileD20, time, fineU_uphi[0]);
		fineMesh.toTecplot(outputFileD21, time, fineU_uphi[0]);
		fineMesh.toTecplot(outputFileD22, time, fineU_uphi[0]);
		
		double[][] ESED = new double[3][fineMesh.getNv()];
		fineMesh.toTecplot(outputFileDDf, time, fineU_uphi[0]);
		fineMesh.toTecplot(outputFileDD0, time, fineU_uphi[0]);
		fineMesh.toTecplot(outputFileDD1, time, fineU_uphi[0]);
		fineMesh.toTecplot(outputFileDD2, time, fineU_uphi[0]);
		double[] myParam = new double[1];
		
		String errorFile = "AxisPED\\RerrorFile_"+N+".dat";
		FileWriter fw = null;
		try {
			fw = new FileWriter(errorFile);
		} catch (IOException e) {
			// TODO: handle exception
		}
		
		double[][] gradient = new double[1][fineMesh.getNv()];
		while(time <= endTime) {
			time += dt;
			myParam[0] = time;
			System.out.println("time: " + time);
			//fs.applyBC_MBN(RHS, Direct.All, 0.0, 1, 2);
			//fsHomo.applyBC_MBN(RHSHomo, Direct.X, 0.0, 1);
			fsHomo.applyBC_MBN(RHSHomo, Direct.Y, 0.0, 1, 2);
			//fsFine.applyBC_MBN(RHSFine, Direct.X, 0.0, 1);
			fsFine.applyBC_MBN(RHSFine, Direct.Y, 0.0, 1, 2);
			//System.out.println(Arrays.toString(RHS));
			
			uHomo = homoSolver.LUSolve(RHSHomo, 0, 0.005);
			uFine = fineSolver.LUSolve(RHSFine, 0, 0.005);
			for(int i = 0; i < fsHomo.GetNdof(); i++) {
			    aHomo[i] = c0 * (uHomo[i] - u0Homo[i]) - c2 * v0Homo[i] - c3 * a0Homo[i];
			    vHomo[i] = v0Homo[i] + c6 * a0Homo[i] + c7 * aHomo[i]; 
			    uvaHomo[i] = c0 * uHomo[i] + c2 * vHomo[i] + c3 * aHomo[i];				
			}
			
			for(int i = 0; i < fsFine.GetNdof(); i++) {
			    aFine[i] = c0 * (uFine[i] - u0Fine[i]) - c2 * v0Fine[i] - c3 * a0Fine[i];
			    vFine[i] = v0Fine[i] + c6 * a0Fine[i] + c7 * aFine[i]; 
			    uvaFine[i] = c0 * uFine[i] + c2 * vFine[i] + c3 * aFine[i];				
			}
			
			MHomo.mul(uvaHomo, RHSHomo);
			MFine.mul(uvaFine, RHSFine);
			for(int i = 0; i < fsHomo.GetNdof(); i++) {
				u0Homo[i] = uHomo[i];
			    v0Homo[i] = vHomo[i];
			    a0Homo[i] = aHomo[i];
			}
			
			for(int i = 0; i < fsFine.GetNdof(); i++) {
				u0Fine[i] = uFine[i];
			    v0Fine[i] = vFine[i];
			    a0Fine[i] = aFine[i];
			}
			
			
			ss.assembleAxisPESOTSSolution(fsFine, fineU0, fineU1, fineU2, fsHomo, uHomo, fsCell, FOCS, SOCS, N, ep);
			fsFine.extractComponent(uFine, fineU_uphi);
			fsFine.extractComponent(fineU0, fineU0_uphi);
			fsFine.extractComponent(fineU1, fineU1_uphi);
			fsFine.extractComponent(fineU2, fineU2_uphi);
			u012[0] = fineU_uphi[0];
			u012[1] = fineU0_uphi[0];
			u012[2] = fineU1_uphi[0];
			u012[3] = fineU2_uphi[0];
			
			phi012[0] = fineU_uphi[1];
			phi012[1] = fineU0_uphi[1];
			phi012[2] = fineU1_uphi[1];
			phi012[3] = fineU2_uphi[1];
			
			fineMesh.toTecplot(outputFile1f, time, fineU_uphi[0]);
			fsFine1.computeGradient(fineU_uphi[0], gradient);
			fineMesh.toTecplot(outputFileD1f, time, gradient[0]);
			
			fineMesh.toTecplot(outputFile10, time, fineU0_uphi[0]);
			fsFine1.computeGradient(fineU0_uphi[0], gradient);
			fineMesh.toTecplot(outputFileD10, time, gradient[0]);
			
			fineMesh.toTecplot(outputFile11, time, fineU1_uphi[0]);
			fsFine1.computeGradient(fineU1_uphi[0], gradient);
			fineMesh.toTecplot(outputFileD11, time, gradient[0]);
			
			fineMesh.toTecplot(outputFile12, time, fineU2_uphi[0]);
			fsFine1.computeGradient(fineU2_uphi[0], gradient);
			fineMesh.toTecplot(outputFileD12, time, gradient[0]);			

			//====================================================
			fineMesh.toTecplot(outputFile2f, time, fineU_uphi[1]);
			fsFine1.computeGradient(fineU_uphi[1], gradient);
			for(int i = 0; i < fineMesh.getNv(); i++) {
				gradient[0][i] = -gradient[0][i];
			}
			fineMesh.toTecplot(outputFileD2f, time, gradient[0]);
			
			fineMesh.toTecplot(outputFile20, time, fineU0_uphi[1]);
			fsFine1.computeGradient(fineU0_uphi[1], gradient);
			for(int i = 0; i < fineMesh.getNv(); i++) {
				gradient[0][i] = -gradient[0][i];
			}
			fineMesh.toTecplot(outputFileD20, time, gradient[0]);
			
			fineMesh.toTecplot(outputFile21, time, fineU1_uphi[1]);
			fsFine1.computeGradient(fineU1_uphi[1], gradient);
			for(int i = 0; i < fineMesh.getNv(); i++) {
				gradient[0][i] = -gradient[0][i];
			}
			fineMesh.toTecplot(outputFileD21, time, gradient[0]);
			
			fineMesh.toTecplot(outputFile22, time, fineU2_uphi[1]);
			fsFine1.computeGradient(fineU2_uphi[1], gradient);
			for(int i = 0; i < fineMesh.getNv(); i++) {
				gradient[0][i] = -gradient[0][i];
			}
			fineMesh.toTecplot(outputFileD22, time, gradient[0]);
			
			fsFine.computeStressAndElectricDisp(matFunc, fineU_uphi, ESED);
			fineMesh.toTecplot(outputFileDDf, time, ESED[0]);

			fsFine.computeStressAndElectricDisp(DHomo, fineU0_uphi, ESED);
			fineMesh.toTecplot(outputFileDD0, time, ESED[0]);

			fsFine.computeStressAndElectricDisp(matFunc, fineU1_uphi, ESED);
			fineMesh.toTecplot(outputFileDD1, time, ESED[0]);

			fsFine.computeStressAndElectricDisp(matFunc, fineU2_uphi, ESED);
			fineMesh.toTecplot(outputFileDD2, time, ESED[0]);
			try {
				double L2Norm = fsFine1.computeL2Norm(fineU_uphi[0]);
				double[] error0 = MVO.add(fineU_uphi[0], -1.0, fineU0_uphi[0]);
				double[] error1 = MVO.add(fineU_uphi[0], -1.0, fineU1_uphi[0]);
				double[] error2 = MVO.add(fineU_uphi[0], -1.0, fineU2_uphi[0]);
				fw.write(time + "\t");
				fw.write(fsFine1.computeL2Norm(error0) / L2Norm + "\t");
				fw.write(fsFine1.computeL2Norm(error1) / L2Norm + "\t");
				fw.write(fsFine1.computeL2Norm(error2) / L2Norm + "\t");
			
				double H1SemiNorm = fsFine1.computeH1SemiNorm(fineU_uphi[0]);
				fw.write(fsFine1.computeH1SemiNorm(error0) / H1SemiNorm + "\t");
				fw.write(fsFine1.computeH1SemiNorm(error1) / H1SemiNorm + "\t");
				fw.write(fsFine1.computeH1SemiNorm(error2) / H1SemiNorm + "\t");
				
				L2Norm = fsFine1.computeL2Norm(fineU_uphi[1]);
		        error0 = MVO.add(fineU_uphi[1], -1.0, fineU0_uphi[1]);
				error1 = MVO.add(fineU_uphi[1], -1.0, fineU1_uphi[1]);
				error2 = MVO.add(fineU_uphi[1], -1.0, fineU2_uphi[1]);

				fw.write(fsFine1.computeL2Norm(error0) / L2Norm + "\t");
				fw.write(fsFine1.computeL2Norm(error1) / L2Norm + "\t");
				fw.write(fsFine1.computeL2Norm(error2) / L2Norm + "\t");
			
				H1SemiNorm = fsFine1.computeH1SemiNorm(fineU_uphi[1]);
				fw.write(fsFine1.computeH1SemiNorm(error0) / H1SemiNorm + "\t");
				fw.write(fsFine1.computeH1SemiNorm(error1) / H1SemiNorm + "\t");
				fw.write(fsFine1.computeH1SemiNorm(error2) / H1SemiNorm + "\n");
				
			} catch (IOException e) {
				// TODO Auto-generated catch block
				e.printStackTrace();
			}

		}
		
		try {
			fw.flush();
			fw.close();
		} catch (IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
		
	}

}
